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Extension of Marcatili’s analytical approach for 
rectangular silicon optical waveguides 
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Abstract —Marcatili’s famous approximate analytical descrip¬ 
tion of light propagating through rectangular dielectric wave¬ 
guides, published in 1969, gives accurate results for low-index- 
contrast waveguides. However, photonic integrated circuit tech¬ 
nology has advanced to high-index-contrast (HIC) waveguides. 
In this paper, we improve Marcatili’s model by adjusting the 
amplitudes of the components of the electromagnetic fields in 
his description. We find that Marcatili’s eigenvalue equation 
for the propagation constant is also valid for HIC waveguides. 
Our improved method shows much better agreement with rig¬ 
orous numerical simulations, in particular for the case of HIC 
waveguides. We also derive explicit expressions for the effective 
group index and the effects of external forces on the propagation 
constant. Furthermore, with our method the phenomenon of 
avoided crossing of modes is observed and studied. 

Index Terms —Optical waveguides. Electromagnetic propaga¬ 
tion, Electromagnetic fields. Integrated optics, Silicon on insula¬ 
tor technology. Optical sensors. 

I. Introduction 

The propagation of light through rectangular dielectric 
optical waveguides cannot be described in closed analytical 
form. Marcatili’s famous approximate analytical approach 
[1] has been used since the 1970’s and is treated in many 
textbooks on optical waveguides theory [2]-[5]. His method is, 
however, derived for waveguides with a low refractive index 
contrast, while nowadays technology has advanced to high- 
index-contrast (HIC) waveguides. Silicon-on-insulator (SOI) 
technology has, for example, become one of the focus plat¬ 
forms for integrated optics over the last decade. The large 
refractive index contrast of the materials allows for small 
device footprint. High-yield mass fabrication is provided using 
CMOS processes from the electronics industry, that have been 
tailored to photonic applications [6]. Behavior of integrated 
optical components, such as ring resonator filters or arrayed 
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waveguide grating (AWG) based multiplexers depend critically 
on the exact knowledge of the propagating modes in the 
waveguide [7], [8]. Although numerical solutions such as the 
circular harmonics method, the film mode matching (FMM) 
method, the variational mode expansion method (VMEM), 
or the finite element method (FEM) are available [9]-[12], 
we believe that an analytical model is useful in order to 
gain insight in the physics of the devices, and also for fast 
explorative simulations of photonic integrated circuitry [13]. 

In this paper, we extend the range of waveguides for which 
Marcatili’s approximate approach can be applied, in particular 
to high-index-contrast waveguides. Similar to Marcatili, we 
use an ansatz for the form of the modal fields that is based 
on separation of variables in the waveguide core. The large 
index contrast causes, in Marcatili’s original approach, a severe 
mismatch of the electromagnetic fields just inside and just 
outside the core of the waveguide. We show that Marcatili’s 
eigenvalue equation for the propagation speed of the light 
through waveguides is, in fact, more general and we obtain 
improved modal electromagnetic fields for the same eigenvalue 
equation, which have a much lower mismatch. An analytical 
description is presented, and is compared with the fundamental 
mismatch of this ansatz, which is found by means of an 
optimization algorithm. 

Next to this, explicit equations are derived for the effective 
group index and for the linearized infiuence of external effects 
on the effective index of the modes. As an example, we 
analytically calculate the infiuence of temperature on the 
effective index of the modes in the waveguide. Also, results 
are presented on photonic evanescent field sensors, where the 
refractive index of the medium in the vicinity of the waveguide 
is probed with the evanescent tail of the waveguide mode [14], 
[15]. 

Throughout this paper, we test the analysis with the first 
three modes in a typical SOI waveguide with a guiding layer 
height of 300 nm, with infrared light that has a free-space 
wavelength Aq = 1550 nm. These guides consist of a thin 
monocrystalline silicon layer (n = 3.476) on top of a thick 
silicon dioxide (BOX) layer (n = 1.444) [16]. The infiuence 
of the silicon substrate below the BOX layer is neglected. 

In the next section, we present our extension of Marcatili’s 
approach. In Sec. Ill, we apply the eigenvalue equation for the 
effective index and derive explicit equations for the effective 
group index and for the effects of external forces. In Sec. IV, 
we compare our theory with rigorous simulations of typical 
SOI waveguides. Section V concludes this paper. 
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11. Theory 

In this section, we describe an approximate model of how 
light travels through rectangular high-index-contrast (HIC) 
dielectric waveguides. In Sec. II-A, we make an ansatz for 
modes in a rectangular waveguide, and briefly investigate the 
implications of this ansatz at the interfaces of the waveguide 
core. In Sec. II-B, we describe Marcatili’s choice of parameters 
for the ansatz and Sec. II-C presents our improved method. In 
Sec. II-D, we present a quantification of the error that arises 
from the discontinuities of the electromagnetic field at the 
interfaces, and we propose an optimization method based on 
this quantification. In Sec. II-E, we use this quantification to 
discuss the different methods. 



A. Rectangular waveguides 

In regular lossless dielectric waveguides, the core of the 
waveguide has a higher refractive index (ni) than the sur¬ 
rounding media (n 2 -n^), as depicted in Fig 1. The refractive 
indices for the outer quadrants, i.e. in the corner regions, are 
not specified because we neglect these regions. We consider a 
monochromatic wave with angular frequency cc, propagating in 
the waveguide direction (z) with a propagation constant (3. For 
a two-dimensional refractive index profile n(x, y), solutions of 
Maxwell’s equations can be found in the form 

S{x,y,z,t) =Re{E{x,y)ex.p[i{ujt- /3z)]}, (1) 


and a similar description of H, with S and H being the 
electric and magnetic field, respectively. The free-space prop¬ 
agation constant is ko = ujjc, where c is the speed of light in 
vacuum. When the propagation constant P is larger than the 
propagation constant that is allowed in the regions outside the 
waveguide core, due to spectral cut-off (i.e. p > koUj, with 
j = 2,.., 5), the light is confined in the core of the waveguide. 
The lateral confinement of such guided waves dictates the light 
to exist in the form of certain modes, or “standing waves”. 

Using Maxwell’s equations, it is now possible to describe 
the full electromagnetic fields in terms of the longitudinal field 
components [2], i.e. in region j 
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where yo and eo are the permeability and the permittivity of 
vacuum, and Kj is defined by 


Ej = - / 32 , ( 6 ) 

All components, and hence in particular the longitudinal 
components Ez and Hz, satisfy the reduced wave equation 


Fig. 1. Schematics of a rectangular waveguide with a width of 600 nm and 
a height of 300 nm, used to guide light with a free-space wavelength of 1550 
nm. (a) 3-D sketch of the waveguide, also displaying the tangential electric 
field components, (b) Cross-section of the waveguide at 2 ; = 0. Regions 1 
to 5 are indicated. The shaded corner regions are neglected in this analysis. 
The color plot represents Ex, the dominant electric field component, of the 
fundamental TM-like mode, (c) Sketch of the cross-section of the waveguide 
with the coordinate frame for a typical TM-like mode. Shape of Ex in the x 
direction (solid line) and in the y-direction (dashed line). Refractive indices 
ni = 3.476, 712 = 1.444, 773 = 714 = 775 = 1. (d) Sketch of the cross- 
section of the waveguide with the coordinate frame for a typical TE-like mode. 
The waveguide geometry is rotated such that Ex is tangential to the “upper” 
surface of the waveguide. Refractive indices 771 = 3.476, 774 = 1.444, 
r72 = 773 = 775 = 1. The mode profiles in this figure are calculated using 
the amplitude optimization method. 


(here given for Ez) [2]: 
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Furthermore, when Fqs. (2)-(5) are satisfied, the electric and 
magnetic fields satisfy V • {eon‘jS) = 0 and V • = 0, 

respectively. 

In this section, we adopt a description of the behavior 
of the light in a rectangular waveguide that is based on 
separation of spatial variables (x, y) in the core region, similar 
to Marcatili’s ansatz. We neglect the effect of the corners 
based on the observation that the field is small in those 
areas. The modal field then consists of standing waves in the 
core of the waveguide and an exponentially decaying field 
outside the core (see Fig 1). We will show that the form we 
propose here can not provide an exact solution of Maxwell’s 
equations, and thus provides an approximate description of the 
physics. The proposed solution obeys Maxwell’s equations in 
regions 1-5, so that all errors that arise from the approximation 
show themselves at the interfaces between the waveguide 
core and its cladding. At these interfaces, continuity of the 
electromagnetic field components tangential to these interfaces 
is required (referred to as the electromagnetic boundary condi¬ 
tions), but we allow for discontinuities that are small compared 
to the field strength in the core of the waveguide. When all 
tangential components are continuous, the normal components 
automatically obey Maxwell’s equations. 

In our analysis, we propose a form where the electric field 
is predominantly polarized in the x-direction. From symmetry, 
the field could as well be chosen predominantly polarized in 
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the y-direction as there is no further discrimination between 
the X- and the y-direction. Sketches (a)-(c) in Fig. 1 show 
the fundamental TM-like mode in a waveguide with a width 
of 600 nm and a height of 300 nm. Sketch (d) shows 
the fundamental TE-like mode in this waveguide. Instead of 
rotating the coordinate-frame, the waveguide itself was rotated 
such that the “top” surface of the waveguide is tangential to 
Ex- 

We make the ansatz on the form of the longitudinal com¬ 
ponents of the modal electromagnetic field in region 1 : 

=Ai sm[k:o{x + 0] cos[fcj;(y + ??)], (8) 

=A 2 cos[k:o{x + C)] sm[ky{y + tj)]. (9) 


Then the transversal electromagnetic field components inside 
region 1 follow from Eqs. (2)-(5) as 
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In region 2, we set 


(13) 


Ez =A3exp[y2{x + d/2)] cos[ky{y + y)], (14) 

Hz =A 4 exp[ 72 (a:; + d/2)] sm[ky{y + y)], (15) 

in region 3 

Ez =A5exp[-'y3{x - d/2)]cos[ky{y + y)], (16) 
Hz =Aeexp[-'y3{x - d/2)]sm[ky{y + y)], (17) 

in region 4 

Ez =Ajsm[kx{x + ^)]exp['y4{y + b/2)], (18) 

Hz =As cos[A:a;(a; + C)] exp[74(y + 6/2)], (19) 


and in region 5 

Ez =Ag sm[kx{x + C)] exp[- 75 (y - 6/2)], (20) 

Hz =Aio cos[kx{x + ^)] exp[- 75 ( 2 / - 6 / 2 )]. ( 21 ) 

Here, amplitudes Ai - Aio, spatial frequencies kx, ky, spatial 
shifts ^, 77 , and exponential decay strengths 7 j > 0 , j = 2,..,5 
are still to be determined. The transversal components are 
calculated from Eqs (2)-(5). From the wave equation (7) in 
region 1 , we arrive at 

Kf = nlkl -I3^ = kl+kl, ( 22 ) 

or, with positive /3, 


/3 = ky- 


The effective (refractive) index of the mode in the waveguide is 
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The sine and cosine dependency on x and y, of the transver¬ 
sal electric (Ez) and magnetic (Hz) field components, was 
chosen such that the dominant electric field component (Ex) 
is described by a cosine function in both x and y direction. 
For most common waveguides 722 ~ 723 and 724 ^ 725 , from 
which follows that spatial shifts ^ and y are small. From 
the observation that the field of the fundamental mode of 
a waveguide has its highest energy density in the center of 
the guide, we expect that the field components with a cosine 
behavior in both x- and y-directions (Ex and Hy) carry the 
majority the field’s energy, next the components with a sine 
and a cosine dependence in x or y (Ez and Hz), and finally 
the least energy is expected in the field components that have 
a sinusoidal profile in both x- and y-direction (Ey and Hx). 
These components have antinodes (i.e. high energy density) 
in the corners of the waveguide core. 

This description has Ai - Aio, 77 , kx and ky as free 
parameters. There are 4 interfaces, with 4 tangential field 
components to be matched per interface, adding up to 16 
equations from the electromagnetic boundary conditions. The 
tangential field components are depicted in Fig. 1(a). Field 
amplitude is used as normalization factor, so we end up 
with an overdetermined system of only 13 free parameters for 
16 equations. 

In summary, we proposed an ansatz on the form of the 
electromagnetic fields of the modes in a rectangular dielectric 
waveguide, such that Maxwell’s equations are obeyed in 
all regions 1-5. This ansatz has 13 free parameters. From 
continuity of the tangential electromagnetic field components 
at the interfaces, 16 boundary condition equations follow. In 
the remainder of Sec. II-A, the requirements that follow from 
continuity at either the surface normal (Sec. II-Al) or parallel 
(Sec. II-A2) to the dominant electric field component, are 
given. 

1) Obeying the electromagnetic boundary conditions at the 
interfaces normal to the dominant electric field component: 
In this section, we derive the requirements that follow from 
continuity of the fields at the 1-2 and 1-3 interfaces. The 
dominant electric field component. Ex , is orthogonal to these 
interfaces, so an infinitely wide, in the y-direction (h 00 ), 
rectangle will describe a TM mode in a slab waveguide, see 
Fig. 1(c). From all eight electromagnetic boundary conditions 
at these interfaces, we find 


4 (xeonjky ^ 

I3k. 
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As =Ai sm[kxiC - d/2)], 

(29) 

A4 =A2 cos[kx{^ - d/2)], 

(30) 

As =Ai sm[kxi^ + d/2)], 

(31) 


(23) 
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Ae =A2 cos[k:c(^ + d/2)\, (32) 


together with 


2 

tan[fca;(^ - (i/2)] = (33) 

77,2 i^x 
2 

iim[k^{i + d/2)\ = (34) 

kx 

Equations (29)-(32) follow directly from the continuity of Ez 
and Hz . The continuity of Ey and Hy is most easily verified by 
substituting Eqs. (28)-(34) into the remaining electromagnetic 
boundary conditions. With these field amplitudes 742 -^ 6 , the 
magnetic field component Hx is zero in regions 1, 2 and 3, 
as follows from Eq. (4). 

The last two equations, (33) and (34), can be recognized as 
the eigenvalue equations for a TM mode in a slab waveguide 
[2], [4]. These eigenvalue equations thus do not only hold for 
a slab solution where d/dy = 0 and Hy, Ex and Ez are the 
non-zero field components, but also for our ansatz where there 
is a variation in the y-direction. We eliminate ^ from these two 
equations, and arrive at the functional 


G{kx,ko,ni,n2,n3,d) = tdin[kxd]- 


nffca;(n|72 + n|73) 


nlnlkl - nf 7273 


= 0 . 

(35) 

2 ) Obeying the electromagnetic boundary conditions at the 
interface parallel to the dominant electric field component: 
When we obey all eight boundary conditions at the 1-4 and 
1-5 interface, to which the dominant electric field component. 
Ex, is parallel, we find 


together with 


I3ky 
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Aj =Ai cos[ky{r] - 6 / 2 )], 
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^8 =A 2 sm[ky{ri - 6 / 2 )], 
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Ag =Ai cos[ky{ri + 6/2)], 
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II 

(41) 

tdin[ky{r] b/2)] 

II 

Cn 

(42) 


Eliminating 77 from latter two equations gives 
F{ky,ko,ni,n 4 ,n 5 ,b) = tan[fcj^ 6 ] - 

ky — 7475 


Equations (37)-(40) follow from the continuity of Ez and 
Hz. Continuity of Ex and Hx is checked by substituting 
Eqs. (36)-(42) into the four boundary conditions corresponding 
to these field components at the two interfaces. It follows 
from Eq. (3) that with these field amplitudes ^ 42 , Aj - ^ 410 , 
the electric field component Ey is zero in regions 1, 4 and 
5. Equations (41) and (42) are identical to the eigenvalue 
equations for a TE mode in a slab waveguide. 


3) Conclusion: The eigenvalue equations that follow from 
the analysis of slab waveguides, in which invariance of the 
field in one direction is assumed, are in fact more general. 
Identical equations follow from the ansatz for Ez and Hz, 
i.e. Eqs. (8)-(21), and from imposing the boundary conditions 
on the horizontal or vertical interfaces of the rectangular 
waveguide. 

Obeying all electromagnetic boundary conditions at the 
interfaces normal to Ex (Sec. II-Al) demand a different 
amplitude coefficient of the magnetic field in the core {A 2 ) 
then the conditions that follow from the interfaces parallel to 
Ex (Sec. II-A2). Therefore the ansatz has no solutions that 
excactly obeys the electromagnetic boundary conditions at all 
interfaces simultaneously. The next sections are devoted to 
different possibilities for choosing the free parameters such 
that a low mismatch of the fields at the boundaries is achieved. 

4) Normalization: Throughout this work, Ai is normalized 
such that the power fiux through waveguide regions 1-5 equals 
unity, i.e. 


P = 



II 

regions 1-5 


EyH*) dxdy > = 1, 


(44) 


where the integral runs over the regions 1-5. 


B. MarcatiWs approach 

Marcatili has developed a widely used analytical approach 
for low-index-contrast waveguides [1]. Eor propagating modes 
in these guides, k^rij « (3 because modes are not guided 
otherwise, so kxlkoUj and kyjkorij are much smaller then 
unity. Therefore those quantities are neglected in second order. 
This is often referred to as “far from cutoff”, while, “close to 
core material spectral cutoff” would be more appropriate. 

In Marcatili’s work, Hx is set to zero in all regions. With this 
requirement, all electromagnetic boundary conditions at the 1 - 

2 and 1-3 interface can be satisfied. This was shown in Sec. 
II-Al, where the requirement that Hx = 0 in regions 1, 2 and 

3 followed from the boundary conditions at these interfaces. 
At the 1-4 and 1-5 interfaces approximations are necessary 
since not all boundary conditions can be satisfied at these 
interfaces. Hx is set to zero and continuity follows trivially. 
In the approximate matching of Ex across these interfaces, 
quantities on the order of {kx/kt^rijY are neglected, which 
is reasonable for low-index-contrast waveguides. Eor these 
guides, y/JEHz is larger than -^/edEz, and is matched across 
the horizontal interfaces, while Ez is not matched. Erom the 
requirements above, it is found that the eigenvalue equation 
of this waveguide is identical to the slab eigenvalue equations 
(33), (34), (41), and (42). We will refer to this approach as 
MarcatiWs Hx = 0 method. 

Although neglecting terms on the order of {kx/konj)‘^ is 
valid for low-index-contrast waveguides, this quantity is even 
larger than unity outside the core region of high-index-contrast 
waveguides. This approximation introduces a large mismatch 
in the continuity of Ex, which is the dominant electric field 
component. 
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Similarly, another approximate solution is obtained by 
setting Ey = 0 [1], [3]. With this demand, all boundary 
conditions at the interface parallel to can be satisfied, and 
mismatches occur at the 1-2 an 1-3 interfaces. Analogue to the 
approach above, Ey is trivially matched and Ez is matched, 
but Hz is not. Hy, which is the dominant magnetic field 
component, is matched while neglecting terms on the order of 
{kx/kQUjY. Although the field amplitudes are different to the 
case where = 0 in all regions, the eigenvalue equations of 
this waveguide description are identical and thus also given by 
Eqs. (33), (34), (41), and (42). We will refer to this approach 
as MarcatiWs Ey = ^ method. 


C. Improved method 

In this section, we present two improvements of Marcatili’s 
methods that give a better description of high-index-contrast 
waveguides. Two choices for the matching of the boundary 
conditions are presented, one where the fields at the interfaces 
normal to the dominant electric field component, E^, are 
continuous, and one where the fields at the interfaces parallel 
to Ex are continuous. In Sec. II-E, we show that the latter is 
more accurate for the cases considered in this paper. 

We argue that the dominant boundary conditions for deter¬ 
mining kx and ^ are at the 1-2 and 1-3 interfaces, and therefore 
Eqs. (33) and (34) are used to determine kx and Similarly, 
continuity at the 1-4 and 1-5 interfaces is used to determine 
ky and rj using Eq. (41) and (42). This is supported by the 
argument that the two different approximations of the modal 
electromagnetic field as presented by Marcatili both yield these 
eigenvalue equations. 

1) Improved Hx ~ 0 method: In comparison with Mar¬ 
catili’s approach, we remove the discontinuity of the dominant 
electric field at the cost of the weak magnetic field component, 
Hx, being not continuous across the 1-4 and 1-5 interfaces. 
In this approach, we demand continuity of the tangential 
electromagnetic fields at the 1-2 and 1-3 interfaces, which 
are normal to Ex. This determines the electromagnetic fields 
in regions 1, 2 and 3 by amplitudes A 2 -Aq, analogue to 
MarcatiWs Hx = 0 method. At the 1-4 and 1-5 interfaces, 
we choose to match Ex and Hz perfectly. The first is the 
dominant electric field component. The latter is chosen in favor 
of matching Ez because an infinitely high (d ^ oo) waveguide 
is identical to a slab waveguide in which Ez = These 
requirements determine the electromagnetic field in regions 4 
and 5, in which Hx is not necessary zero. We will refer to 
this method as the Improved Hx ^ ^ method. 

This gives the slab eigenvalue equations (35) and (43), and 
field amplitude Eqs. (28)-(32), together with 

Aj =Ai cos[fcj^(7? - b/2)], (45) 

As =A 2 sin[fcj^(r? - b/2)], (46) 

AQ=Ai(^ + ^A!h—J!E^cos[ky{r] + b/2)], (47) 

Aio =A 2 s'm[ky{ri + b/2)]. (48) 


2) Improved Ey ^ D method: This improvement removes 
the discontinuity of the dominant magnetic field component 
Hy, that was present in MarcatiWs Ey = ^ method. We match 
all tangential field components at the 1-4 and 1-5 interface, 
which are parallel to the dominant electric field component. 
Erom Sec. II-A2 follows that = 0 in regions 1, 4, and 5. 
At the 1-2 and 1-3 interfaces, Hy is matched because it is 
the dominant magnetic field component, and Ez is matched 
in favor of Hz because an infinitely wide waveguide (h oc) 
is identical to a slab waveguide in which Hz = 0. 

We find the slab eigenvalue equations (35) and (43), and 


field amplitude Eqs. (36)-(40), together with 

^3 =Ai sin[fca;(^ - d/ 2 )], (49) 

A 4 =A 2 (^1 + cos[kx{^ - d/2)], (50) 

As =^isin[fc^(e + d/2)], (51) 

Ae = A 2 (^1 + MMzLI) ^ cos[fc, + d/2)]. (52) 


D. Least-discontinuity optimization of the ansatz parameters 

In Sec. II-A, we presented an ansatz on the form of the 
electromagnetic field for modes in a rectangular waveguide. 
This ansatz was chosen such that Maxwell’s equations are 
satisfied in regions 1-5, so that all errors manifest themselves 
at the four interfaces of the waveguide core. This error, 
which is the discontinuity of the tangential electromagnetic 
field components at these interfaces, is referred to as the 
mismatch. The measure we adopt to quantify this mismatch, 
or error, is the average energy density that is associated with 
these discontinuities. In Sec. IV, we show that this intuitive 
quantity excellently agrees with rigorous numerical results. 
This analysis is performed in a cross-section of the waveguide 
at 2 ; = 0, at time t = 0, as further longitudinal and temporal 
behavior follows trivially from Eq. (1). We define: 

C/mm = I ^(n++n-)2.|i>x (£;+-£;-) fdl (53) 

+ ^^|i> X {H+ -H-fd\. 

The four interfaces of the waveguide are simultaneously de¬ 
scribed by the integral. The line integral runs along the circum¬ 
ference of the waveguide in the (x,y)-plane, and I = 2(6-fd) is 
the length of this circumference. and E~ are the electric 
fields just outside and inside the waveguide core region 1, so 
that {E~^—E~) represents the discontinuity of this field, and 0 
is a unit vector orthogonal to the waveguide surface. The cross 
product of 0 with the discontinuity in the field just selects the 
tangential components. and n~ are the refractive indices 
just outside and inside the waveguide. At the interface, an 
average refractive index (n+-Fn“)/2 is assumed to calculate 
the energy density of the electric field components. 

Although f/mm can be intuitively interpreted as an energy 
density, we cannot attach a rigorous physical meaning to this 
quantity. The mismatch in the fields only occurs at interfaces, 
which have no physical volume. Therefore the energy density 
cannot be integrated over volume in order to obtain a total 
energy. 
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In this section, we use this quantity to propose a new (semi-) 
analytical method. In the next section, we use this quantity to 
investigate the methods of the two previous sections (II-B and 
II-C) without resorting to numerical simulations. 

In what we call the full optimization method, the mismatch 
in the fields at the boundaries, as given by Eq. (53), is 
minimized using an unconstrained nonlinear optimization as 
implemented in Matlab [17]. The initial parameters are cal¬ 
culated from the improved method as presented in Sec. II-C2. 

The full optimization method may be compared with the 
variational mode expansion method (VMEM) [11], which is 
also based on a separation of variables. However, the methods 
differ substantially because the modal field ansatz differs and 
also the object function used in the variational principle of 
VMEM differs from the object function that is minimized 
in our approach. The VMEM is applicable to more general 
structures. 

It is interesting to know how accurate and ky, and 
thus the effective refractive index, are calculated using the 
analytical methods. The effective indices that follow from the 
eigenvalue equations on the one hand, and from the least- 
mismatch optimization on the other hand, are presented in Eig. 
2. It can be seen that the difference between the two methods 
is small. The influence of this difference in propagation speed 
on the mismatch of the fields at the interfaces is investigated by 
optimizing the field amplitudes A 2 - Aio with fixed kx, ky, 
and T] as calculated from the analytical eigenvalue equations. 
As Eq. (53) is quadratic in the amplitudes A 2 - Aio, this 
minimum can be found analytically. This method is referred 
to as the amplitudes optimization method. As can be seen 
in Eig. 3, the mismatch of the method with the fitted /c’s 
(kx and ky) and method with the analytical calculated k's is 
almost identical. Therefore we conclude that, with the ansatz 
for the fields as described in Sec. II-A, and the error of the 
model described by the field mismatch, Eq. (53), the values 
of kx and ky are very accurately calculated from analytical 
slab eigenvalue equations, for the typical SOI waveguides as 
considered in this paper. 

E. Discussion of the different methods 

In Eig. 3, all six different methods are compared. The energy 
associated with mismatch of the electromagnetic fields at the 
boundaries of the core of the waveguide, f/mm, is plotted for 
three types of waveguide modes, f/mm can be interpreted as 
the energy density of the error at the interfaces. In order to 
get a feeling for the magnitude of f/mm, the average energy 
density in the core region of the waveguide is also plotted in 
this figure (black dashed line). 

Eor all geometries, the improved methods are better than 
Marcatili’s original methods. Considering the fundamental 
modes in a 400 nm wide waveguide, and the 1st order TE-like 
mode in a 650 nm wide waveguide, we see that the mismatch 
in Marcatili’s Ey = 0 method is 1.5 to 4 times larger than 
the mismatch in the improved Ey ^ ^ method. In Marcatili’s 
method, the power fiux through regions 2, 3, 4, or 5, i.e. in 
the evanescent tail of the mode, is sometimes negative (or in 
backwards propagating direction). This is an indication of the 



waveguide width [nm] 


Fig. 2. Effective refractive indices calculated using four different methods. 
Plot (a) presents the first 3 modes in the waveguide core (TEO, TMO, TEl). In 
comparison with conventional notation, e.g. TEqo, we dropped a zero, as in 
our waveguide geometries all higher-order modes have higher-order standing 
waves only in the direction of the width of the waveguide. The numerically 
calculated effective index of the TMl mode is included for completeness. 
Plots (b)-(d) show the difference in the effective index as calculated by the 
analytical method with respect to the numerical method. 

inaccuracy of the method, and this behavior was not observed 
for the other methods. 

Eor the waveguide geometries as considered in this paper, 
the improved Ey ^ method has smaller errors in the 
continuity of the fields at the boundaries than the improved 
Hx ^ ^ method. However, this is geometry dependent, as 
is expected from the trend in the TM-like mode. Eor TM- 
like modes in waveguides wider than 1765 nm (h > Gd), the 
improved Hx ~ 0 method has a lower mismatch than the 
improved Ey ^ method. This behavior of the asymptotes is 
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400 500 600 700 800 900 1000 1100 


waveguide width [nm] 

— €> - - Marcatili Hx=0 Amplitudes optimization 

— X- - - Improved Hx~0 . Full optimization 

O Mareatili Ey=0 —^— Fitting to numerieal profile 

-X - Improved Ey~0 (used in Fig. 4 only) 

Fig. 3. Mismatch of the tangential field components at the interfaces of the 
waveguide core, Umm, see Eq. (53). Typical, 300 nm high, SOI waveguide. 
Plots (a)-(c) represent the TEO, TMO, and TEl modes in the waveguide. In 
each plot, 6 different methods to calculate the free parameter in the analytical 
description of the waveguide mode profile are presented. The black dashed 
line indicates the average energy density in the core (calculated using the full 
optimization method). 


expected from theory of slab waveguides. In the limit 6 —> oo, 
the waveguide is a slab waveguide with interfaces 1-2 and 
1-3, so the improved ^ ^ method is expected to give 
better results, because the fields described by this methods 
are continuous across these interfaces. In the limit d —> oo, 
only the 1-4 and 1-5 interfaces play a role, so we expect the 
improved Ey ^ method to give the most accurate results. 
For the cases considered in this plot, average energy density 
associated with the mismatch (f/mm) of the improved Ey ^ 
method is below 3% of the average energy density in the core 
of the waveguide. 

In conclusion, the effective index of the mode is accurately 
calculated by the eigenvalue equations. The amplitudes A 2 to 
A 10 are calculated more accurately with the improved methods 
than with Marcatili’s original methods. A better quantitative 
evaluation can be made by a comparison of the methods 
with rigorous numerical calculated modal profiles, which is 
presented in Sec. IV. 

III. Novel applications of the eigenvalue equation 

OF THE PROPAGATION CONSTANT 

In this section, we derive explicit expressions for the ef¬ 
fective group index and for the infiuence of a change in the 
waveguide. 


A. Modal dispersion 

The modal dispersion can be described by the (effective 
refractive) group index: 

dn^ii 8(3 
dko' 


— TT-eff A — 


dXo 


(54) 


The group index is often used to design photonic integrated 
circuits. For example to calculate of the free spectral range 
(FSR) of ring resonators. From Eq. (23), we find 


dp 

dko 


drii 


= — konf -h - k 




(55) 


dkx 

dko dko 

The 1st and 2nd term on the right-hand-side of this equation 
are specified by the material refractive indices. The refrac¬ 
tive indices nj{ko) may depend on frequency and thus on 
/co = cc/c. The 3rd term is calculated from Eq. (35). Although 
kx is only given as implicitly, dkx/dko can be calculated 
explicitly. The total derivative of the left-hand-side of Eq. 
(35) with respect to ko, dGjdko, equals zero for solutions 
of G = 0. The height d does not depend on frequency. So we 


get 






dG 

_ 9G 

dG 

1 _ 

dk^ dG 

dni 

dko 

dko 

dkx 

dko drii 

dko 

or, 







dkx _ 

0:1 Q: 

, dG 
' dni 

dni 1 
dko ' 

dG 

dn2 


8G dn2 

dn 2 dko 


dG dns 

dns dko 


dko 


dG 

dk^ 


8G 8ns 

dns dko ’ 

(56) 


Similarly, the 4th term of the right-hand-side of Eq. (55) is 
calculated from Eq. (43) as 


dky 

dko 


dF 

dko 


dF dm 
dni dko 


dF dn4 
dn4 dko 


dF driB 
dn 5 dko 


dF 

dk^, 


(57) 


The partial derivatives in Eqs. (56) and (57) are straightforward 
to calculate, as shown in appendix A. 
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B. Influence of a linear external effect 

Analogue to the derivation of the effective group index, it is 
also possible to calculate the linearized influence of an external 
effect on the propagation constant of the mode. For example, 
the influence of a temperature change in the waveguide or the 
influence of a change in the refractive index of the surrounding 
medium can be calculated. We describe the external effect by a 
parameter %. In this theory, the waveguide properties rij, h and 
d are assumed to be known, as well as the flrst-order influence 
of the external effect on these properties, i.e. drijldx^ ddjdx 
and dhjdx- Taking the derivative of Eq. (23) with respect to 
X gives 


dx 


>9ni 


dK 


dx 


dx 


i I 


(58) 


As in the previous section, dkxjdx is calculated by taking the 
derivative of Eq. (35), G = 0, and solving dkxjdx' 


dkx 

dx 


^ I I dG dd 

dni dx ' dn2 dx ' dn^ dx ' dd dx 
aG 

dka, 


(59) 


And for dkyjdx, we And from Eq. (43) that 


dky 

dx 


d^drn I d^dr]^ , _^dn5 , dF db 
dni dx ' dn4^ dx dn^ dx ' db dx 
dF 

dky 


(60) 


Equations (59) and (60) are given explicitly in Appendix B. 

1) Influence of temperature on SOI waveguides: In the case 
of a temperature change, drij /dx is given by the thermo-optic 
effect. The change in cross-section of the waveguide, ddjdx 
and dhjdx are described by linear thermal expansion. The 
thermo-optic coefficients of silicon and silicon dioxide, as well 
as the linear thermal expansion coefficients, are reported in 
Refs. [16] and [18], respectively. 

2) Evanescent field sensor: In evanescent held sensors, 

liquids with varying refractive indices flow over the SOI 
waveguides. In this case, x describes the liquid refractive 
index, so that for TM-like modes ns = = x 

and dnsjdx — dn/^jdx — dn^jdx — while the other 
waveguide properties are constant. For TE-like modes, 722 = 
ns = ns = X and dn 2 ldx = dnsjdx = dnsjdx = 1 (see 
Fig. 1). 


IV. Comparison with rigorous numerical methods 

In this section, we compare the analytical results of the 
previous section with rigorous numerical results. To And the 
transverse modes of the waveguide, we use the film mode 
matching (FMM) method as implemented in the FimmWave 
software package [10], [12]. In this method the cross-section 
of the ridge waveguide is split in vertical slices, and 1- 
dimensional modes are computed analytically for each slice. 
The 2-D modes are now found by simultaneously solving 
the modal amplitude factors of the 1-D modes in all slices 
such that they construct a held obeying Maxwell’s equations. 
The area of the numerical simulation extends 2 pm from the 
waveguide, and 200 1-D modes are used per slice. The width 
of the waveguide is varied in steps of 10 nm, and in steps 
of 1 nm between 820 nm and 840 nm. For verification, we 


compared the FMM method with the finite element method 
(FEM) that is also implemented in the FimmWave software, 
using ^210 gridpoints in both the x- and the y-direction. For 
all datapoints presented in this work, the difference (between 
FMM and FEM) in effective index is below 10“^ and the 
relative energy in the difference held is below 10“^. 

The measure that is used to compare two electromagnetic 
fields is the relative energy in the difference held, i.e. 

ff (n^eolE^-E^l^+/iiolH^-H^l^)dxdy 

. regions 1-5 

ff in^eo\EN\^ + lio\H^\^)dxdy 

regions 1-5 

( 61 ) 

where and are the analytically and rigorous numer¬ 
ically calculated fields, respectively. This integral runs over 
all regions that are described by the analytical solution. The 
integral is calculated analytically with the assumption that the 
numerical held is piecewise-constant in space. 


A. Improved analytical methods 

The effective refractive indices are calculated using four 
different methods, including the FMM method, and are pre¬ 
sented in Fig. 2. For completeness, results of the effective 
index method (EIM) [2] are also included. The analytical 
eigenvalue equations, described in Sec. II-A, give, for most 
geometries, the most accurate effective index. The method 
where the effective index is calculated using least-mismatch 
fitting in the boundary conditions as presented in Sec. II-D 
is less accurate for TM-like modes. The error in the EIM 
becomes significantly large for strong confined modes. The 
discontinuities in plots (c) and (d) at a width of 833 nm are 
discussed in Sec. IV-B. The other small “wiggles” in i\\t full 
optimization method are due to numerical instabilities. It can 
be seen that the error in the effective index calculated using the 
analytical eigenvalue equation remains below 2% for all cases. 
The error is, for all methods, lower for the weaker-confined 
modes. 

Figure 4 presents the energy of the difference held between 
the analytical and the numerical FMM method. First of all, it is 
clear that something special happens for the fundamental TM- 
like (TMO) and the 1st TE-like (TEl) modes in the waveguide 
with a width of ^830 nm. We will explain this in Sec. IV-B. 
The very small “wiggles” in the fundamental TE-like mode 
between waveguide widths of 820 nm and 840 nm are due 
to our numerical discretization of the electromagnetic fields, 
which not exactly match the refractive index distribution for 
features below 10 nm. 

The fundamental limit of the ansatz is found by optimizing 
all free parameters such that the energy in the difference held 
compared to the numerical result (the quantity plotted in this 
figure) is minimal. 

The method that gives the least difference with the rigorous 
numerical result is the method where kx, ky, and p follow 
from the slab eigenvalue equations, while the held amplitudes 
A 2 -A 10 are optimized by minimizing the mismatch, f/mm, in 
Eq. (53) {amplitude optimization method). Eor this method, the 
energy in the difference field is below 2% for all waveguide 
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waveguide width [nm] 

Fig. 4. Energy in the electromagnetic difference field between the analytically 
calculated fields and the rigorous numerically calculated fields, normalized to 
the energy in the numerically calculated field. The fields are only considered 
in regions 1-5. Legend and layout of the figure is identical to Fig. 3. 


modes considered here (except for the crossing at 830 nm). 
Also, this method is very close to the fundamental limit of the 
ansatz. 

The fact that the eigenvalue equations work so well, justifies 
the assumption that the boundary conditions at the 1-2 and 1-3 
interfaces can be best used to calculate kx and the boundary 
conditions on the 1-4 and 1-5 interfaces can be best used 
to calculate ky. Deviating from these /c-values to minimize 
the total mismatch in the boundary conditions only slightly 
reduces this mismatch and increases the error in the /c-values. 
Still, the quantification of the mismatch as done in Sec. II-D 
gives surprisingly accurate results. The difference between 
the fundamental limit of the ansatz and both the amplitude 


optimization method and full optimization method is very 
small, especially considering that there is a fundamental error 
in the ansatz, and that the error of the fields only manifests 
itself at the interfaces of the waveguide core, which is difficult 
to interpret from a physical perspective. 

The improved approaches are, for all cases, more accurate 
then Marcatili’s original approaches. For these methods, the 
full field-profiles are easily and intuitively calculated. The 
fields of the SOI waveguides considered in this work are 
accurately described by the improved Ey ^ ^ method, with 
less than 3% energy in the difference field (except for the 
crossing at 830 nm). 

The improved Ey ^ ^ method is more accurate than the 
improved Hx ^ ^ method for (h < 3d), i.e. for TE-like modes 
and for TM-like modes in waveguides that are not very wide. 
From the mismatch in the boundary conditions, the same trend 
was expected and this crossing was expected at (b ^ 6d). 
Also, the difference between the improved Ey ^ ^ method, 
the improved Hx ^ method and the optimization methods 
follows the same trend as expected from Sec. II-D. 

In conclusion, the improved Marcatili’s method describes 
the effective indices and the field distributions of typical 
SOI waveguides sufficiently accurately for many applications, 
except when the effective indices of two modes are similar. In 
the next section, the modes found from our approximation are 
used to explain what happens in these regions. When compared 
to rigorous numerical solutions, the error in the effective index 
was below 2% and the relative error in terms of the energy 
of the fields was less than 3%. In addition, the optimization 
method gives even more accurate results for the fields. 

B. The avoided crossing phenomenon 

Figures 2 and 4 suggest that something interesting happens 
at the apparent crossing of the effective indices of the TMO- 
like and the TEl-like modes when the width of the guide is 
changed. A detailed inspection of the waveguides with widths 
around 833 nm is presented in Fig. 5. In plot (a), it can be seen 
that the numerically computed effective indices of the 2nd and 
3rd mode in the waveguide (counted from high to low effective 
index) actually do not cross each other, but show a behavior 
that is known in quantum mechanics as avoided crossing [19]. 
We investigate the modes that where found numerically, E 2 
and , in terms of the analytically computed approximate 
modes of the waveguide. 

We will verify that the actual modes E 2 and E^ can 
in good approximation be written as a superposition of the 
TMO-like, F^tmo, and TEl-like, F^tei modes that where cal¬ 
culated using the approximate theory {amplitude optimization 
method), i.e. 

E^ ^ aEjuo + ^^TEi, (62) 

for some real a and b, and i = 2 or 3. The phase of mode Ei 
is chosen such that coefficient b is positive. The coefficient a 
of the TMO-like mode can be either positive or negative. The 
approximate calculated modes F^tmo and F^tei, are in good 
approximation orthonormal such that normalization of Ei in 
the norm of Eq. (44) implies b = Vl — o?. 
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Fig. 5. Investigation of the avoided crossing of effective indices of two 
modes. Plot (a): numerically calculated effective indices of the 2nd and 3rd 
mode, zoom-in of Fig. 2(a). Plot (b): power in the TMO-like mode when the 
fields of the modes in plot (a) are written as a superposition of a TMO-like 
and a TMl-like mode. The curves in plot (a) are color-coded accordingly. 
Plot (c): Relative energy in the electromagnetic difference field between the 
superposition and the rigorous numerically calculated fields, to be compared 
with Fig 4 (b,c). 


The coefficient a of the TMO-like mode is optimized such 
that the difference measured using Eq. (61) between the left- 
and right-hand sides of Eq. (62) is minimum. The result is 
plotted in Eig. 5(b), where it can be seen that mode 2 looks 
like a TMO-mode at the left of the crossing, while it looks 
like a TEl-like mode on the right-hand-side of the crossing, 
whereas close to the crossing the modes are an equal mixture 
of E^tmo and E^tei- In Fig- 5(c), it can be seen that the error 
between the superposition Ei and the rigorous numerically 
calculated field close to the apparent crossing is small 
and similar to the error that was found away from the crossing 
(see Eig. 4). Therefore we may indeed conclude that the 
field around the crossing can be written as a superposition of 
modes of the types that are present away from the crossing. 
Eigure 6 presents the electric fields E^tmo, ^tei, E 2 and 
E^ for a 833 nm wide by 300 nm high waveguide, where 
a‘^ ^ ^ 0.5. 

Using this observation, we will derive a qualitative de¬ 
scription of this avoided crossing. The exact description of 
propagation of light through waveguides that are invariant in 



y y y 


Fig. 6. Electric fields of different modes in a 833 nm wide by 300 nm 
high waveguide. The dashed lines separate the different regions (see Fig. 1). 
Modes Ejmo, and E^ are plotted from top to bottom. Color 

indicates the field strength, white regions have a low field strength. Regions 
with positive field strength are indicated with a plus (+) sign and red. Regions 
with a negative field strength are indicated with a minus sign (-) and blue. 


the z-direction can be formulated as an eigenvalue problem 
with the propagation constant as eigenvalue. In fact, by 
eliminating the longitudinal components and from 
Maxwell’s equations, one obtains the following eigenvalue 
problem for the transverse components only 


O 


Ex{x,y)\ 

Ey(x,y) 

Hx{x,y) 

\Hy{x,y)J 


Pi 


Ex{x,y)\ 

Ey{x,y) 

Hx{x,y) 

\Hy{x,y)J 


(63) 


with the propagation constant /3i as eigenvalue and with O 
a second order partial differential operator with respect to 
transverse variables x and y. We consider forward propagating 
waves with positive f3i. The operator O is not symmetric, 
however it can be shown that solutions of Eq. (63) are 
orthogonal with respect to the bilinear form derived from the 
power flux [3] 


00 

{i\j) = jj [Ei X H*) ■ z dxdy = 6ij (64) 

— OO 


where the inner product between two solutions i and j is 
defined, 6ij is the Kronecker delta function, and the eigenfields 
have been normalized. We adopt a bra-ket notation and write 
Eq. (63) as 

d\i) = m ( 65 ) 


We now apply the aforementioned observation that in good 
approximation the electromagnetic fields of the modes in the 
waveguide, also around the crossing, can be written as a 
superposition of the approximate fields, hence 


Ei ^ aJE^TMO + ^^TEi, or \i) ^ a\a) + 6|5), (66) 
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where \a) and \b) represent the TMO-like and TEl-like modes 
in the waveguide, respectively, while \i) represents the exact 
solutions of Eq. (65). We only consider the 2nd and 3rd 
approximate solutions here. As will become clear, only modes 
with similar propagation constants have to be taken into 
account around the crossing, whereas the other modes are 
already accurately calculated by the approximate methods 
presented in this paper. Substituting Eq. ( 66 ) in Eq. (65) and 
taking the inner product with (a| gives 

a(a|Oa) + b{a\Ob) ^ Pi {a{a\a) + b{a\b )). (67) 


If we also take the inner product of Eq. (65) with {b\ we arrive 
at the ( 2 x 2 )-system: 


/(a|Oa) (a| 06 )A /a\ 

[{b\6a) {b\6b)) \b) 


A 


f{a\a) 


{a\b)\ 

mj 



or, 



(69) 


where M = 

J_ / {b\b){a\da) — {a\b){b\da) {b\b){a\db) — {a\b){b\db) \ 

D y — {b\a) {a\Oa) + {a\a){b\da) —{b\a){a\db) + {a\a){b\db) J 


and 

D = {a\a){b\b) — {a\b){b\a). 


The modes that we found in our approximate analysis are 
almost orthonormal, so (a|a) and {b\b) are approximately unity 
and (a16) and {b\a) are approximately zero. Away from the 
crossing, we found that the approximate solutions \a) and \b) 
obey relation Eq. (65) so that (a|Oa) « Pa, {b\Ob) ^ P^, 
while (a 106) and ( 6 |Oa) are small. This allows us to write 
the relation Eq. (69) as 


f Pa ^ab ^ 

V ha Pb-\- hj \bJ 



(70) 


where Sa, Sab, and S^a are quantities that are much smaller 
than the P's. This system has the eigenvalues [19] 


P2,3 


PL + Pb , V{/3'a-/3ir+iSatSta 
2 2 


with corresponding eigenvectors ^ 2,3 (not normalized) 


(71) 


2Sab _ 

~Pa + (P'a — P'P)^ + ^SabSba 

where P'^ = Pa ^ Sa and + 65 . The two propagation 

constants are closest when p'^ = p'^ but are always separated 
by a minimum distance 4:^/6abSba, so that they never intersect. 
Eor small Sa^Si,^ Sab-, Sba ^ \Pa — Pb\, we find the eigenvector 
for Pa > Pb to be V 2 ~ (1,0) and ^3 « (0,1)- The upper 
propagation constant, P 2 , has a TMO-like mode in this limit, 
while the lower propagation constant, P^, has a TEl-like mode. 
Eor Pb > pa we find V 2 ~ (0,1) and ^3 (1, 0), so that the 

upper propagation constant now has a TEl-like mode while 
the lower propagation constant has a TMO-like mode. An 
interesting case occurs when p'^ = Pl^ and Sab = Sba- Then the 
normalized eigenvectors of this system are • ( 1 , 1 ) 

and t ;3 = ^-(l,— 1 ), i.e. they are either an addition or 
subtraction of the eigenvectors far from the crossing. 



This simple description of the avoided crossing agrees with 
the observations of the numerically computed modal profiles 
as presented in Eigs. 5 and 6 . 

C. Novel applications of the eigenvalue equation of the prop¬ 
agation constant 

We applied the analytical eigenvalue equations for the 
effective index of the modes in rectangular waveguide to 
explicitly calculate of the effective group index (Sec. III-A), 
as well as the linearized influence of an external effect (Sec. 
III-B). In this section, we compare those analytical results with 
rigorous numerical simulations. Eor numerical calculation of 
the linearized infiuence of an external effect on the effective 
index of the modes in the waveguide, we use the following 
procedure: first, we calculated the effective refractive index 
(rieff) for different values of the external effect (x) that are 
around the starting value (xo), then, a linear fit is used to 
obtain the dn^f^/dx- 

1) Effective group index: The method presented in Sec. 
III-A is investigated in this section. Equations (54), (77) and 
(78) are used to calculate the effective group index of the 
modes in the waveguides. The silicon dispersion is taken into 
account, with dni/dko = 3.147 • 10^ m“^ [16]. These results 
are compared with the group index that was numerically calcu¬ 
lated using the EMM method as implemented in EimmWave. 
Results are presented in Eig. 7 where can be seen that the 
error remains below 4%. 

2) Temperature-induced effective index change: To rigor¬ 
ously compute the linearized infiuence of temperature varia¬ 
tions on the effective index of the waveguide, the temperature 
in the simulation is varied from 20 °C to 28 °C in steps of 2 
°C. The infiuence dn^a/dT is compared with the analytically 
calculated infiuence in Eig. 8 . Two effects are taken into 
account: the thermo-optic effect using dnxjdx — 1-83 • 10“"^ 

and the thermal expansion using dbjdx — m/^C and 
ddldx = Oid m/°C, with linear thermal expansion coefficient 
(Y = 2.6 • 10“^ [16], [18]. The effects are investigated 

simultaneous (net effect) and the effect of the expansion only 
is also shown separately. This separate investigation does, of 
course, not simulate a physical measurement but it gives a 
good test case for the method. Erom Eig. 8 , it can be concluded 
that the effects are well estimated with a relative error below 
7%. The error is much smaller for modes that are less confined, 
and the effect of the geometrical change is better estimated 
than the effect of the silicon index change. 

3) Cladding-index induced effective index change: In Sec. 
III-B2, we described the linearized infiuence of a change in 
refractive index of the cladding of the waveguide. This is, for 
example, used in evanescent field sensors, where a liquid fiows 
over the waveguide surface. Eor example, aqueous solutions 
with different concentrations of sucrose or sodium chlorine are 
used in Refs. [14] and [15], respectively. The rigorous numer¬ 
ical results are obtained by changing the cladding refractive 
index from 1.310 (water) to 1.320 in steps of 0.002. Erom this, 
5neff/9nciad is extracted. Results are shown in Eig. 9, where 
it can be seen that relative errors are up to 23%. These errors 
rapidly decrease for wider waveguides, i.e. for less-confined 
modes. 
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waveguide width [nm] 


Fig. 7. Effective group indices. Plots (a), (c), (e) present the numerical 
result, as well as the result of the extension of Marcatili’s method (eigenvalue 
equation). Plots (b), (d), (f) show the difference between the two methods. 


Si index change & expansion Expansion only 

index change, dmddT [10 "^/K] index change, dn^ddl [10 ^/K] 



waveguide width [nm] waveguide width [nm] 


Eig. 8. Shift in effective refractive index due to temperature variation. (Left 
plots) Shift in effective refractive index due to change in silicon refractive 
index (thermo-optic effect) as well as expansion of the waveguide (thermal 
expansion) is depicted. (Right plots) Only the geometry of the waveguide 
is varied (thermal expansion), with the silicon refractive being constant. The 
result of the extension of Marcatili’s method, and the result of the numerical 
EMM method are presented. 


index change, dnQwdnc\?i& [-] index change, dnQwdnc\a& [-] 




waveguide width [nm] 


Eig. 9. Shift in effective index due to change in cladding refractive index. 
The unit is refractive index change per cladding material index change [-/-]. 
The initial cladding is water (n = 1.31). 
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V. Conclusion 

We derived a model for high-index-contrast rectangular 
dielectric waveguides, based on an ansatz that describes the 
fields in the core of the waveguide as standing waves. As in 
Marcatili’s work, we also arrive at eigenvalue equations for the 
spatial frequencies of the standing waves, i.e. and ky, that 
are identical to the eigenvalue equations of slab waveguides. In 
order to obtain kx, we use the limit that the waveguide extends 
to infinity in the y-direction, and order to obtain ky, we use the 
limit that the waveguide extends to infinity in the x-direction. 
The effective refractive index of the mode in the waveguide, 
which directly follows from kx and ky, agrees excellently with 
rigorous numerical mode solvers (relative error < 2%). In 
comparison with Marcatili’s original work, our novel choice 
of electromagnetic field amplitudes in the ansatz removed the 
discontinuity of the dominant electromagnetic field compo¬ 
nents, which was severe for high-index-contrast waveguides. 
This improvement led to better agreement with numerical 
simulations, with a relative energy of the difference field below 
3%, except when the effective indices of two modes in the 
waveguide are similar. We quantified the error that arises from 
the discontinuity of the fields in terms of energy density. With 
this quantification, we accurately predicted the trend in the 
error of the different analytical methods. In addition, this 
quantification allowed us to optimize the amplitudes in the 
ansatz such that minimal mismatch is achieved. Our amplitude 
optimization method accurately describes both TE-like as well 
as TM-like modes in the rectangular waveguides (except when 
two modes have similar effective indices), and the error of this 
method with respect to the numerical results is very close to 
the fundamental error in the ansatz. 

Next to this, we derived explicit expressions for the effective 
group index of the waveguide, taking dispersion into account. 
The error in comparison with rigorous numerical methods 
was below 4%. Also, explicit expressions are derived for the 
linearized influence of an external effect. We predicted the 
infiuence of temperature, as well as the influence of a change 
of the refractive index of the cladding of the waveguide. 

We applied our method to interpret the results of a rigorous 
numerical mode solver, and showed that the modes in the 
waveguide show avoided crossing behavior when the effective 
indices of a TE-like and a TM-like mode are similar. A quick 
analysis of waveguides can be performed with this method, 
while the results of this method can be used as starting point 
for rigorous mode solvers. Marcatili’s method has long been 
used for educational purposes, as is indicated by its frequent 
occurrence in textbooks. With our work, the development of 
intuitive analytical methods now follows the technological 
developments to high-index-contrast waveguides. 

Appendix A 

Explicit equation eor modal dispersion 

In Sec. Ill-A, partial derivatives dkx/dko and dkx/dko 
were given by Eqs. (56) and (57). In this section, we give 
dkx/dko and dky/dko for the case when only the core ma¬ 
terial is dispersive, i.e. ni{ko) and where the other refractive 
indices do not depend on the frequency, thus also not on ko. 
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And find that Eqs. (56) and (57) are explicitly given by Eqs. 
(77) and (78), which are shown on the next page. 


Appendix B 

Explicit equation for external effect 

In this appendix, we give the explicit form of Eqs. (59) and 
(60). Let: 
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Then Eqs. (59) and (60) are identical to Eqs. (81) and (82), 
see next page, where a 2 -OL^ are defined in Eqs. (73) - (76). 
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